use "$data/NLS_hh_data.dta", clear

***** Table 3: Horse race between migration and covid symptoms

* Collapse pre-data hunger in April (in case surveyed twice)
recode series (1=0)
keep if month==11 // | month==10
collapse village migr_any highMigV symptoms wgt_c (max) hungry , by(hhid month series)
	// take max of pre-hunger

* Get April--April difference in hunger
bys hhid: egen wgt = max(wgt_c)
drop wgt_c
reshape wide hungry, i(hhid month) j(series)
gen byte h_diff = hungry2 - hungry0

* Run regressions (no upazila FEs - want to keep cross-upazila COVID variation)
label var migr_any "Migrant Household"
label var highMigV "High Migration Village"
label var symptoms "COVID Symptoms"

eststo clear
eststo: reg h_diff migr_any [aw=wgt], vce(cluster hhid)
eststo: reg h_diff migr_any symptoms [aw=wgt], vce(cluster hhid)
eststo: reg h_diff highMigV [aw=wgt], vce(cluster hhid)
eststo: reg h_diff highMigV symptoms [aw=wgt], vce(cluster hhid)
eststo: reg h_diff migr_any highMigV symptoms [aw=wgt], vce(cluster hhid)

local col1 " & \multicolumn{5}{c}{Change in Food Insecurity April 2019--April 2020} \\"
local col2 " & (1) & (2) & (3) & (4) & (5) \\ \midrule "	
esttab using "$tables\Table_3.tex", keep(migr_any highMigV symptoms) order(migr_any highMigV symptoms) ///
	replace se tex label  mlabels(none) nonotes nonumbers nostar b(3) /// starlevels(* 0.1 ** 0.05 *** 0.01)
	stats(r2 N, fmt(2 0) labels("R-Squared" "Observations")) ///
	posthead(`col1' `col2')


***** Monthly event study graphs (Figures 5A and 6B)
use "$data/NLS_hh_data.dta", clear

recode series (1=0)

* Split by calendar month, NLS/COVID survey, and HH migration status
reg hungry month#series#migr_any [aw=wgt_c], vce(cluster hhid)
predict bM, xb
predict seM, stdp

* Split by calendar month NLS/COVID survey, migration status, and village migration status
reg hungry month#series#migr_any#highMigV [aw=wgt_c], vce(cluster hhid)
predict bV, xb
predict seV, stdp

collapse bM seM bV seV, by(month series migr_any highMigV) fast

* Alternative: collapse (mean) b_h=hungry (sem) se_h=hungry [aw=wgt_c], by(.)
* This gets the same coefficients as the regression, sem is close to clustered regression standard errors

***** Make graphs for Figs. 5 and 6

* Duplicate first/last month for graphing
expand 2 if month==1, gen(newmig)
replace month = 13 if newmig
drop newmig 

* Offset and jiggle the dots a little
replace month = month-0.5 
replace month = month+.09 if series==2 & migr_any==1

foreach i in M V {
	gen hi`i' = b`i' + 1.96*se`i'
	gen lo`i' = b`i' - 1.96*se`i'
}

** Figure 5
local covdate 9
tw 	(scatteri -.025 `covdate' .41 `covdate', recast(line) lcolor(gs0) lpattern(shortdash)) ///
	(rarea loM hiM month if series==0 & migr_any==0 & highMigV==0, color(gs0%15) lwidth(0)) ///
	(rarea loM hiM month if series==0 & migr_any==1 & highMigV==0, color(gs0%15) lwidth(0)) ///
	(line bM month if series==0 & migr_any==0 & highMigV==0, color(gs0))  ///
	(line bM month if series==0 & migr_any==1 & highMigV==0, color(red*.8))  ///
	(rcap loM hiM month if series==2 & migr_any==0 & highMigV==0, color(gs0) msize(small) lwidth(thin)) ///
	(rcap loM hiM month if series==2 & migr_any==1 & highMigV==0, color(red*.8) msize(small) lwidth(thin)) ///
	(scatter bM month if series==2 & migr_any==0 & highMigV==0, mcolor(gs0) msymbol(O) )  ///
	(scatter bM month if series==2 & migr_any==1 & highMigV==0, mcolor(red*.8) msymbol(T))  ///
	, graphregion(color(white)) ytitle("Share of Households that are Food Insecure") ///
	xlabel(1(1)12, angle(45) format(%tdm) )   ///
	 ylabel(0(0.1)0.4, angle(0)) xtitle("")   ///
	 ttext(.42 `covdate' "COVID-19")  plotregion(margin(0)) ///
	 xlabel(1 "Jul" 2 "Aug" 3 "Sep" 4 "Oct" 5 "Nov" 6 "Dec" 7 "Jan" 8 "Feb" 9 "Mar" 10 "Apr" 11 "May" 12 "Jun", angle(45)) ///
	 legend(order (9 "Has Migrant, 2020" 8 "No Migrant, 2020" 5 "Has Migrant, Typical Year" 4 "No Migrant, Typical Year"))  

graph export "$figures\Figure_5A.png", replace

* Jiggle the dots some more
replace month = month+.09 if series==2 & migr_any==1
replace month = month-.09 if series==2 & highMigV==0
//replace month = month+.05 if series==2 & migr_any==1 & highMigV==0
//replace month = month+.1 if series==2 & migr_any==1 & highMigV==1

** Figure 6
local covdate 9
tw 	(connected bV month if series==. & migr_any==0 & highMigV==0, lpattern(solid) msymbol(O) color(gs0)) ///
	(connected bV month if series==. & migr_any==0 & highMigV==1, lpattern(dash) msymbol(Oh) color(gs0)) ///
	(connected bV month if series==. & migr_any==1 & highMigV==0, lpattern(solid) msymbol(T) color(red*.8)) ///
	(connected bV month if series==. & migr_any==1 & highMigV==1, lpattern(dash) msymbol(Th) color(red*.8)) ///
	(scatteri -.05 `covdate' .57 `covdate', recast(line) lcolor(gs0) lpattern(shortdash)) ///
	(rarea loV hiV month if series==0 & migr_any==0 & highMigV==0, color(gs0%15) lwidth(0)) ///
	(rarea loV hiV month if series==0 & migr_any==0 & highMigV==1, color(gs0%15) lwidth(0)) ///
	(rarea loV hiV month if series==0 & migr_any==1 & highMigV==0, color(gs0%15) lwidth(0)) ///
	(rarea loV hiV month if series==0 & migr_any==1 & highMigV==1, color(gs0%15) lwidth(0)) ///
	(line bV month if series==0 & migr_any==0 & highMigV==0, color(gs0))  ///
	(line bV month if series==0 & migr_any==1 & highMigV==0, color(red*.8) lpattern(solid))  ///
	(line bV month if series==0 & migr_any==0 & highMigV==1, color(gs0) lpattern(dash))  ///
	(line bV month if series==0 & migr_any==1 & highMigV==1, color(red*.8) lpattern(dash))  ///
	(rcap loV hiV month if series==2 & migr_any==0 & highMigV==0, color(gs0) msize(small) lwidth(thin)) ///
	(rcap loV hiV month if series==2 & migr_any==1 & highMigV==0, color(red*.8) msize(small) lwidth(thin)) ///
	(rcap loV hiV month if series==2 & migr_any==0 & highMigV==1, color(gs0) msize(small) lwidth(thin)) ///
	(rcap loV hiV month if series==2 & migr_any==1 & highMigV==1, color(red*.8) msize(small) lwidth(thin)) ///
	(scatter bV month if series==2 & migr_any==0 & highMigV==0, mcolor(gs0) msymbol(O))  ///
	(scatter bV month if series==2 & migr_any==1 & highMigV==0, mcolor(red*.8) msymbol(T))  ///
	(scatter bV month if series==2 & migr_any==0 & highMigV==1, mcolor(gs0) msymbol(Oh))  ///
	(scatter bV month if series==2 & migr_any==1 & highMigV==1, mcolor(red*.8) msymbol(Th))  ///
	, graphregion(color(white)) ytitle("Share of Households that are Food Insecure") ///
	xlabel(1(1)12, angle(45) format(%tdm) )   ///
	 ylabel(, angle(0)) xtitle("")  ///
	 ttext(.585 `covdate' "COVID-19")  plotregion(margin(0)) ///
	 xlabel(1 "Jul" 2 "Aug" 3 "Sep" 4 "Oct" 5 "Nov" 6 "Dec" 7 "Jan" 8 "Feb" 9 "Mar" 10 "Apr" 11 "May" 12 "Jun", angle(45)) ///
	 legend(order(3 "Has Migrant, Low Village" 1 "No Migrant, Low Village" 4 "Has Migrant, High Village" 2 "No Migrant, High Village"))  

graph export "$figures/Figure_6B.png", replace
